%% NOTE DO NOT USE THIS. THIS GIVES WRONG ANSWER AS COMPARED TO THE FORTRAN PROGRAM BY SCHERLIESS
% vdriftmodel(9.333, -76, [268;172.828 ]) GIVES A DRIFT  19.8347 M/S
% where as the actual value, according to the fortran program is 18.07 !!
% manoj nair May , 2012
% vdrift-model_jgr0499.0 -- translated by f2c (version 20100827).
% System generated locals */

%       ************************************************************ */
%       ************************************************************ */
%       SUBROUTINE CALCULATES EQUATORIAL VERTICAL DRIFT AS DESCRIBED */
%       IN SCHERLIESS AND FEJER, JGR, 104, 6829-6842, 1999 */
%       ************************************************************ */
%       INPUT:   XT: SOLAR LOCAL TIME */
%                XL: GEOGRAPHIC LONGITUDE (+ EAST) */

%             PARAM: 2-DIM ARRAY (DOY,F10.7CM) */
%                    DOY     :Day of Year has to run from 1 to 365 (366) */
%                    F10.7cm : F10.7cm solar flux */

%       OUTPUT:   Y: EQUATORIAL VERTICAL DRIFT */
%       ************************************************************ */
% Parameter adjustments */
% Table of constant values */

function[y] = vdriftmodel(xt, xl, param)

% Initialized data */

coeff = [-10.80592, -9.63722, -11.52666, -.05716, ...
    -.06288, .03564, -5.80962, -7.86988, -8.50888, -.05194, -.05798, ...
    -.00138, 2.09876, -19.99896, -5.11393, -.0537, -.06585, .03171, ...
    -10.22653, -3.62499, -14.85924, -.04023, -.0119, -.09656, ...
    -4.8518, -26.26264, -6.20501, -.05342, -.05174, .02419, ...
    -13.98936, -18.10416, -9.30503, -.01969, -.03132, -.01984, ...
    -18.36633, -24.44898, -16.69001, .02033, -.03414, -.02062, ...
    -20.27621, -16.95623, -36.58234, .01445, -.02044, -.08297, ...
    1.4445, 5.53004, 4.55166, -.02356, -.04267, .05023, 5.50589, ...
    7.05381, 1.94387, -.03147, -.03548, .01166, 3.24165, 10.05002, ...
    4.26218, -.03419, -.02651, .07456, 7.02218, .06708, -11.31012, ...
    -.03252, -.01021, -.09008, -3.47588, -2.82534, -4.17668, -.03719, ...
    -.01519, .06507, -4.02607, -11.19563, -10.52923, -.00592, ...
    -.01286, -.00477, -11.47478, -9.57758, -10.36887, .04555, ...
    -.02249, .00528, -14.19283, 7.86422, -8.76821, .05758, -.02398, ...
    -.04075, 14.5889, 36.63322, 27.57497, .01358, -.02316, .04723, ...
    12.53122, 29.38367, 21.40356, -7.1e-4, -.00553, .01484, 18.64421, ...
    26.27327, 18.32704, .00578, .03349, .11249, 4.53014, 6.15099, ...
    7.41935, -.0286, -.00395, -.08394, 14.29422, 9.77569, 2.85689, ...
    -.00107, .04263, .10739, 7.17246, 4.40242, -1.00794, 8.9e-4, ...
    .01436, .00626, 7.75487, 5.01928, 4.36908, .03952, -.00614, ...
    .03039, 10.25556, 8.82631, 24.21745, .05492, -.02968, .00177, ...
    21.86648, 24.03218, 39.82008, .0049, -.01281, -.01715, 19.18547, ...
    23.97403, 34.44242, .01978, .01564, -.02434, 26.30614, 14.22662, ...
    31.16844, .06495, .1959, .05631, 21.09354, 25.56253, 29.91629, ...
    -.04397, -.08079, -.07903, 28.30202, 16.80567, 38.63945, .05864, ...
    .16407, .07622, 22.68528, 25.91119, 40.45979, -.03185, -.01039, ...
    -.01206, 31.98703, 24.46271, 38.13028, -.08738, -.0028, .01322, ...
    46.67387, 16.80171, 22.7719, -.13643, -.05277, -.01982, 13.87476, ...
    20.52521, 5.22899, .00485, -.04357, .0997, 21.46928, 13.55871, ...
    10.23772, -.04457, .01307, .06589, 16.18181, 16.0296, 9.28661, ...
    -.01225, .14623, -.0157, 18.16289, -1.5823, 14.54986, -.00375, ...
    -8.7e-4, .04991, 10.00292, 11.82653, .44417, -.00768, .1594, ...
    -.01775, 12.15362, 5.65843, -1.94855, -.00689, .03851, .04851, ...
    -1.25167, 9.05439, .74164, .01065, .03153, .02433, -15.46799, ...
    18.23132, 27.4532, .00899, -1.7e-4, .03385, 2.70396, -.87077, ...
    6.11476, -8.1e-4, .05167, -.08932, 3.21321, -1.06622, 5.43623, ...
    .01942, .05449, -.03084, 17.79267, -3.44694, 7.10702, .04734, ...
    -.00945, .11516, .46435, 6.78467, 4.27231, -.02122, .10922, ...
    -.03331, 15.31708, 1.70927, 7.99584, .07462, .07515, .08934, ...
    4.19893, 6.01231, 8.04861, .04023, .14767, -.04308, 9.97541, ...
    5.99412, 5.93588, .06611, .12144, -.02124, 13.02837, 10.2995, ...
    -4.862, .04521, .10715, -.05465, 5.26779, 7.09019, 1.76617, ...
    .09339, .22256, .09222, 9.1781, 5.27558, 5.45022, .14749, .11616, ...
    .10418, 9.26391, 4.19982, 12.6625, .11334, .02532, .18919, ...
    13.18695, 6.06564, 11.87835, .26347, .02858, .14801, 10.08476, ...
    6.14899, 17.62618, .09331, .08832, .28208, 10.75302, 7.09244, ...
    13.90643, .09556, .16652, .22751, 6.70338, 11.97698, 18.51413, ...
    .15873, .18936, .15705, 5.68102, 23.81606, 20.65174, .1993, ...
    .15645, .08151, 29.61644, 5.49433, 48.90934, .7071, .40791, ...
    .26325, 17.11994, 19.6538, 44.8881, .4551, .41689, .22398, 8.457, ...
    34.54442, 27.25364, .40867, .37223, .22374, -2.30305, 32.0066, ...
    47.75799, .02178, .43626, .30187, 8.98134, 33.0182, 33.09674, ...
    .33703, .33242, .41156, 14.27619, 20.70858, 50.10005, .30115, ...
    .3257, .45061, 14.44685, 16.14272, 45.40065, .37552, .31419, ...
    .30129, 6.19718, 18.89559, 28.24927, .08864, .41627, .19993, ...
    7.70847, -2.36281, -21.41381, .13766, .05113, -.11631, -9.07236, ...
    3.76797, -20.49962, .03343, .0863, .00188, -8.58113, 5.06009, ...
    -6.23262, .04967, .03334, .24214, -27.85742, 8.34615, -27.72532, ...
    -.08935, .15905, -.03655, 2.77234, .14626, -4.01786, .22338, ...
    -.04478, .1865, 5.61364, -3.82235, -16.72282, .26456, -.03119, ...
    -.08376, 13.35847, -6.11518, -16.50327, .28957, -.01345, -.19223, ...
    -5.3729, -.09562, -27.27889, .00266, .22823, -.35585, -15.29676, ...
    -18.36622, -24.62948, -.31299, -.23832, -.08463, -23.37099, ...
    -13.69954, -26.71177, -.19654, -.18522, -.20679, -26.33762, ...
    -15.96657, -42.51953, -.13575, -.00329, -.28355, -25.4214, ...
    -14.14291, -21.91748, -.2096, -.19176, -.32593, -23.36042, ...
    -23.89895, -46.0527, -.10336, .0303, -.21839, -19.46259, ...
    -21.27918, -32.38143, -.17673, -.15484, -.11226, -19.06169, ...
    -21.1324, -34.01677, -.25497, -.16878, -.11004, -18.39463, ...
    -16.11516, -19.55804, -.19834, -.23271, -.25699, -19.93482, ...
    -17.56433, -18.58818, .06508, -.18075, .02796, -23.64078, ...
    -18.77269, -22.77715, -.02456, -.12238, .02959, -12.44508, ...
    -21.06941, -19.36011, .02746, -.16329, .19792, -26.34187, ...
    -19.78854, -24.06651, -.07299, -.03082, -.03535, -10.71667, ...
    -26.04401, -16.59048, .0285, -.0968, .15143, -18.40481, -23.3777, ...
    -16.3145, -.03989, -.00729, -.01688, -9.68886, -20.59304, ...
    -18.46657, .01092, -.07901, .03422, -.06685, -19.2459, -29.35494, ...
    .12265, -.24792, .05978, -15.32341, -9.0732, -13.76101, -.17018, ...
    -.15122, -.06144, -14.68939, -14.82251, -13.65846, -.11173, ...
    -.1441, -.07133, -18.38628, -18.94631, -19.00893, -.08062, ...
    -.14481, -.12949, -16.15328, -17.40999, -14.08705, -.08485, ...
    -.06896, -.11583, -14.50295, -16.91671, -25.25793, -.06814, ...
    -.13727, -.12213, -10.92188, -14.10852, -24.43877, -.09375, ...
    -.11638, -.09053, -11.64716, -14.9202, -19.99063, -.14792, ...
    -.08681, -.12085, -24.09766, -16.14519, -8.05683, -.24065, ...
    -.05877, -.23726, -25.18396, -15.02034, -15.50531, -.12236, ...
    -.0961, -.00529, -15.27905, -19.36708, -12.94046, -.08571, ...
    -.0956, -.03544, -7.48927, -16.00753, -13.02842, -.07862, -.1011, ...
    -.05807, -13.06383, -27.98698, -18.80004, -.05875, -.03737, ...
    -.11214, -13.6737, -16.44925, -16.12632, -.07228, -.09322, ...
    -.05652, -22.61245, -21.24717, -18.09933, -.05197, -.07477, ...
    -.05235, -27.09189, -21.85181, -20.34676, -.05123, -.05683, ...
    -.07214, -27.09561, -22.76383, -25.41151, -.10272, -.02058, ...
    -.1672];



% Function Body */
index_t = 13;
dim_t = 78;
index_l = 8;
dim_l = 48;
index = 104;
dim = 624;
nfunc = 6;
funct = gfun(param, xl);


%       **********************************
y=0;
%       **********************************
for i=1:index_t
    for il=1:index_l
        kk = index_l*(i-1)+il;
        for j=1:nfunc
            ind = nfunc*(kk-1)+j;
            bspl4=bspl4_time(i,xt)*bspl4_long(il,xl);
            y=y+bspl4*funct(j)*coeff(ind);
        end
    end
end

end % vdrift_model */

% ------------------------------------------------------------------ */
%       ************************************************* */

%       ************************************************* */
function[ret_val] =  bspl4_time(i,  x)
% Initialized data */

t_t = [0., 2.75, 4.75, 5.5, 6.25, 7.25, 10., 14., ...
    17.25, 18., 18.75, 19.75, 21., 24., 26.75, 28.75, 29.5, 30.25, ...
    31.25, 34., 38., 41.25, 42., 42.75, 43.75, 45., 48., 50.75, ...
    52.75, 53.5, 54.25, 55.25, 58., 62., 65.25, 66., 66.75, 67.75, ...
    69., 72];

b = zeros([20,20]);
order = 4;

if(i >= 0)
    if (x < t_t(i-0))
        x=x+24;
    end
end

for j=i:i+order-1
    if(x >= t_t(j) && x < t_t(j+1))
        b(j,1)=1;
    else
        b(j,1)=0;
    end
end

for j=2:order
    for k=i:i+order-j
        b(k,j)=(x-t_t(k))/(t_t(k+j-1)-t_t(k))*b(k,j-1);
        b(k,j)=b(k,j)+(t_t(k+j)-x)/(t_t(k+j)-t_t(k+1))*b(k+1,j-1);
    end
end

ret_val = b(i,order);

end % bspl4_time */

% ------------------------------------------------------------------ */
%       ************************************************* */

%       ************************************************* */
function[ret_val] = bspl4_long(i, x)
% Initialized data */

t_l = [0., 10., 100., 190., 200., 250., 280., 310., ...
    360., 370., 460., 550., 560., 610., 640., 670., 720., 730., 820., ...
    910., 920., 970., 1e3, 1030., 1080];

order = 4;
b = zeros([20,20]);

if(i>=0)
    if (x<t_l(i-0))
        x=x+360;
    end
end
for j=i:i+order-1
    if(x>=t_l(j) && x<t_l(j+1))
        b(j,1)=1;
    else
        b(j,1)=0;
    end
end

for j=2:order
    for k=i:i+order-j
        b(k,j)=(x-t_l(k))/(t_l(k+j-1)-t_l(k))*b(k,j-1);
        b(k,j)=b(k,j)+(t_l(k+j)-x)/(t_l(k+j)-t_l(k+1))*b(k+1,j-1);
    end
end

ret_val=b(i,order);

end % bspl4_long */


function[funct] = gfun (param,  x)


%       *************************************************
flux=param(2);
if(param(2) <= 75)
    flux=75.;
end;
if(param(2) >= 230)
    flux=230.;
end;
cflux=flux;

a=0.;
if((param(1) >= 120) && (param(1) <= 240)) ,a=170.;end
if((param(1) >= 120) && (param(1) <= 240)), sigma=60;end
if((param(1) <= 60) || (param(1) >= 300)) ,a=170.;end
if((param(1) <= 60) || (param(1) >= 300)), sigma=40;end

if((flux <= 95) && (a ~= 0))
    gauss=exp(-0.5*((x-a)^2)/sigma^2);
    cflux=gauss*95.+(1-gauss)*flux;
end
%       *************************************************

%       *************************************************
for i=1:6
    funct(i)=0;
end;
%       *************************************************

%       *************************************************
if((param(1) >= 135) && (param(1) <= 230)), funct(1)=1;end
if((param(1) <= 45) || (param(1) >= 320)), funct(2)=1;end
if((param(1) > 75) && (param(1) < 105)), funct(3)=1;end
if((param(1) > 260) && (param(1) < 290)), funct(3)=1;end
%       *************************************************

if((param(1) >= 45) && (param(1) <= 75))   % W-E
    funct(2)=1.-(param(1)-45.)/30.;
    funct(3)=1-funct(2);
end
if((param(1) >= 105) && (param(1) <= 135))  % E-S
    funct(3)=1.-(param(1)-105.)/30.;
    funct(1)=1-funct(3);
end
if((param(1) >= 230) && (param(1) <= 260))   % S-E
    funct(1)=1.-(param(1)-230.)/30.;
    funct(3)=1-funct(1);
end
if((param(1) >= 290) && (param(1) <= 320))   % E-W
    funct(3)=1.-(param(1)-290.)/30.;
    funct(2)=1-funct(3);
end

%       *************************************************
funct(4)=(cflux-140)*funct(1);
funct(5)=(cflux-140)*funct(2);
funct(6)=(flux-140)*funct(3);

end
% gfun */

